Generic mechanism for generating a liquid-liquid phase transition. 
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Recent experimental results |jj indicate that 
phosphorus, a single-component system, can have 
two liquid phases: a high-density liquid (HDL) 
and a low-density liquid (LDL) phase. A first- 
order transition between two liquids of different 
densities J2]] is consistent with experimental data 
for a variety of materials J^J^], including single- 
component systems such as water [pHf!, silica 
fl9|] and carbon 
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Molecular dynamics simula- 
tions^ qf_very specific models for supercooled wa- 

and supercooled sil- 




11], liquid carbon 
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predict a LDL-HDL critical point, but a 
coherent and general interpretation of the LDL- 
HDL transition is lacking. Here we show that 
the presence of a LDL and a HDL can be directly 
related to an interaction potential with an attrac- 
tive part and two characteristic short-range repul- 
sive distances. This kind of interaction is common 
to other single-component materials in the liquid 



state (in particular liquid metals [14-21,2]), and 
such potentials are often used to decribe systems 
that exhibit a density anomaly [Q. However, our 
results show that the LDL and HDL phases can 
occur in systems with no density anomaly. Our 
results therefore present an experimental chal- 
lenge to uncover a liquid-liquid transition in sys- 
tems like liquid metals, regardless of the presence 
of the density anomaly. 

Several explanations have been developed to under- 
stand the liquid-liquid phase transition. For example, the 
two-liquid models [Q assume that liquids at high pressure 
are a mixture of two liquid phases whose relative con- 
centration depends on external parameters. Other ex- 
planations for the liquid-liquid phase transition assume 
an anisotropic potential |^],[ll-13|. Here we shall see 
that liquid-liquid phase transition phenomena can arise 
solely from an isotropic pair interaction potential with 
two characteristic lengths. 

For molecular liquid phosphorus P4 (as for water), a 
tetrahedral open structure is preferred at low pressures 
P and low temperatures T, while a denser structure is 
favored at high P and high T 0,0,p[. The existence of 
these two structures with different densities suggests a 
pair interaction with two characteristic distances. The 
first distance can be associated with the hard-core exclu- 
sion between two particles and the second distance with 
a weak repulsion {soft-core), which can be overcome at 
large pressure. Here we will use a generic three dimen- 
sional (3D) model composed of particles interacting via 



an isotropic soft-core pair potential. Such isotropic po- 
tentials can be regarded as resulting from an average 
over the angular part of more realistic potentials, and 
are often used as a first approximation to understand 
the qualitative behavior of real systems [fl4}-|22Pl . For 
Ce and Cs, Stell and Hcmmcr proposed a potential with 
nearest-neighbor repulsion and a weak long-range attrac- 
tion p5[ . By means of an exact analysis in ID, they 
found two critical points, with the high-density critical 
point interpreted as a solid-solid transition. Then ana- 
lytic calculations , simulations and exact solution 
in ID [[l8| of the structure factor for a model with a soft- 
core potential were found consistent with experimental 
structure factors for liquid metals such as Bi. The struc- 
ture factor was also the focus of a theoretical study of a 
family of soft-core potentials, for liquid metals, by means 
of mean-spherical approximation | [l9|] ■ More recently, the 
analysis of the solid phase of a model with a soft-core 
potential |^0| was related to the experimental evidence 
of a liquid- liquid critical point in the K2CS metallic al- 
loy pi) ]. Moreover, for simple metals, soft-core potentials 
have been derived by first-principle calculations |l4[ ] and 
computed from experimental data Q . The importance of 
soft-core potentials in the context of supercooled liquids 
was pointed out by Debenedetti and others |2^-|25|, and 
the analysis of experimental data for water gives rise to 
a soft-core potential |2(J . 

The isotropic pair potential considered here (Fig^Ja, 
inset) has a hard-core radius a and a soft-core radius 
b > a. For a < r < b the particles repel each other with 
energy Ur > 0, for b < r < c they attract each other with 
energy — V a < 0. For r > c the interaction is considered 
negligible and is approximated to zero. The potential has 
three free parameters: b/a, c/a and Ur/Ua- 

To select the parameters for the molecular dynam- 
ics (MD) simulations is not easy, and MD is too time 
consuming to study a wide range of parameter values. 
Hence we first perform an integral equation analysis [p7| , 
whose predictions we can calculate very rapidly and effi- 
ciently. In this technique, one derives the phase diagram 
by studying the static pair distribution function - which 
measures the probability of finding a particle at a given 
distance from a reference particle, and thereby quanti- 
fies the correlation between pairs of particles. Mathe- 
matically, this function can be written as the sum of an 
infinite series of many-dimensional integrals, over parti- 
cle coordinates, involving the pair interaction potential. 
Since this exact expression is intractable, approximations 
must be made. Here we use the hypernetted-chain ap- 
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proximation, which consists of neglecting a specific class 
of these integrals, leading to a simplified integral equa- 
tion that can be solved numerically. In the temperature 
- density (T-p) phase diagram, the region where the sim- 
plified equation has no solutions is related |27| to the 
region where the system separates in two fluid phases. 
Thus this technique allows us to estimate the parame- 
ter range where two critical points occur, and hence to 
find useful parameters values for the MD simulations: 
b/a = 2.0, c/a = 2.2 and Ur/Ua = 0.5. We also use in 
the MD calculations several additional parameter sets. 

Specifically we perform MD simulations in 3D at con- 
stant volume V and number of particles N — 490 and 
850. We use periodic boundary conditions, a standard 
collision event list algorithm 24 , and a modified Berend- 
sen method to control T |28| . We find, for each set of pa- 
rameters, the appearance of two critical points (Fig.El). 

A critical point is revealed by the presence of a region, 
in the P-p phase diagram, with negative-slope isotherms. 
In MD simulations this region is related to the coex- 
istence of two phases 0. The (local) maximum and 
minimum along an isotherm correspond to the limits 
of stability of the existence of each single phase (super- 
cooled and superheated phase, respectively). By defini- 
tion, these maxima and minima are points on the spin- 
odal line for that temperature. Since the spinodal line 
has a maximum at a critical point, a way to locate a 
critical point is to find this maximum. In our simula- 
tions (Fig]]]), we find two regions with negatively-sloped 
isotherms and the overall shape of the spinodal line has 
two maxima, showing the presence of two critical points, 
C\ and C%. Using the Maxwell construction in the P- 
V plane Q, we evaluate the coexistence lines of the two 
fluid phases associated with each critical point (Fig]^). 
Considering both the maxima of the spinodal line and 
the maxima of the coexistence regions in the P-p and P- 
T planes, we estimate the low-density critical point C\ at 
Ti = 0.606 ±0.004 U A /k B , Pi = 0.0177 ± 0.0008 U A /a 3 , 
pi = 0.11 ± 0.01 a~ 3 and the high-density critical point 
C 2 at T 2 = 0.665±0.005 U A /k B , P 2 = 0.10±0.01 U A /a 3 , 
P2 = 0.32 ± 0.03 a~ 3 . Critical point C\ is at the end of 
the phase transition line separating the gas phase and 
the LDL phase, while critical point C 2 is at the end of 
the phase transition line separating the gas phase and the 
HDL phase. Their relative positions resemble the phos- 
phorus phase diagram, except that, in the experiments, 
C 2 has not been located jj], but is expected at the end 
of the gas-HDL transition line. 

Our phase diagram (Fig.|^) shows the following fluid 
phases. At high T > T 2 , the only fluid phase is the gas. 
At Xi < T < T 2 , we find - depending on p - the gas 
alone, or the HDL alone (turquoise region), or the HDL 
coexisting with gas (black line in Fig.^Ja and green re- 
gion m Fig]|c). Below Ti, the LDL phase appears alone 
(blue region), or in coexistence with the HDL (orange 
line in Fig.^3 and orange region in Fig]^d), or in coexis- 



tence with the gas (red line in Fig.||b and red region in 
Fig.^d). The point where the gas-LDL coexistence line 
merges with the LDL-HDL coexistence line is the triple 
point. Below the pressure and temperature of the triple 
point, LDL is not stable and separates into gas and HDL. 

For phosphorus the liquid-liquid transition occurs in 
the stable fluid regime 0]. In contrast, for our model, 
it occurs in the metastable fluid regime (see Fig.l). We 
therefore whish to understand how to enhance the sta- 
bility of the critical points with respect to the crystal 
phase. We find that by increasing the attractive well 
width (c — b)/a, both critical temperatures T± and T 2 in- 
crease, and hence both critical points move toward the 
stable fluid phase, analogous to results for attractive po- 
tentials with a single critical point ^,0. For example, 
for attractive well width (c — b)/a = 0.2, both C\ and 
C 2 are metastable with respect to the crystal, while for 
(c — b)/a> 0.7 we find C\ in the stable fluid phase. 

The phase diagram depends sensitively also on the rel- 
ative width of the shoulder b/a and on its relative height 
Ur/Ua- By decreasing b/a or by increasing Ur/Ua, T 2 
decreases and becomes smaller than T±. This means that, 
in these cases, the high-density C 2 occurs below the tem- 
perature of the gas-liquid critical point, i.e. C 2 is in the 
liquid phase and represents a LDL-HDL critical point, as 
in supercooled water f|,|7],[ll| . 

The soft-core potential with the sets of parameter we 
use displays no "density anomaly" (dV/dT)p < 0. This 
result is at first sight surprising since soft-core potentials 
have often been used to explain the density anomaly (see, 
e.g., Refs. p|,|2"4|] ) . To understand this result, we consider 
the entropy S (the degree of disorder in the system) and 
the thermodynamic relation —(dV/dT)p — (dS/dP)T = 
{dS/dV) T {dV/dP) T - Since of necessity (dV/dP) T < 0, 
(dV/dT)p < implies (dS/dV) T < 0, i.e. the den- 
sity anomaly implies that the disorder in the system in- 
creases for decreasing volume. For example, this is the 
case for water. This is consistent with the negative slope 
dP/dT of the crystal-liquid transition line for water, that 
implies AS/AV < for the Clausius-Clapeyron equa- 
tion dP/dT = AS/AV, where AS and AV are the en- 
tropy and volume differences between the two coexisting 
phases. 

For our system, we expect the reverse: (dV/dT)p > 
so (dS/dV)T > 0, consistent with the positive slope of 
the LDL-HDL transition line dP/dT (see Fig]|b). Wc 
confirm our expectation that (dS/dV)T > by explicitly 
calculating S for our system by means of thermodynamic 
integration. 

Our results show that the presence of two critical 
points and the occurrence of the density anomaly are 
not necessarily related, suggesting that one might seek 
experimental evidence of a liquid-liquid phase transition 
in systems with no density anomaly. In particular, a sec- 
ond critical point may also exist in liquid metals that 
can be described by soft-core potentials. Thus the class 
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of experimental systems displaying a second critical point 
may be broader than previously hypothesized. 
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FIG. 1. Pressure-density isotherms, crystallization line 
and spinodal line from the MD simulations for the isotropic 
pair potential in 3D. (a), Inset: The pair potential energy 
U(r) as a function of the distance r between two particles, (a) 
Several isotherms for (bottom to top) ksT /Ua = 0.57, 0.59, 
0.61, 0.63, 0.65, 0.67 (ks is the Boltzmann constant). Dia- 
monds represent data points and lines are guides for the eyes. 
The solid line connecting local maxima and minima along the 
isotherms represents the spinodal line. The two maxima of 
the spinodal line (squares) represent the two critical points 
Ci and Ci. To determine the crystallization line (gray line) - 
below which the fluid is metastable with respect to the crystal 
- we place a crystal seed, prepared at very low T, in contact 
with the fluid, and check, for each (T,p), if the seed grows or 
melts after 10 6 MD steps. The spontaneous formation (nucle- 
ation) of the crystal is observed, within our simulation times 
(« 10 5 MD steps), only for p > 0.27 a" 3 . We use the structure 
factor S(Q) - the Fourier transform of the density-density cor- 
relation function for wave vectors Q - to determine when the 
nucleation occurs. Indeed, at the onset of nucleation, S(Q) 
develops large peaks at finite Q (Q = 12a" 1 and Q — 6a ). 
For each p, we quench the system from a high-T configuration. 
After a transient time for the fluid equilibration, we compute 
P(T,p), averaging over 10 s - 10 6 configurations generated 
from up to 12 independent quenches, making sure that the 
calculations are done before nucleation takes place, (b) En- 
larged view of the region around the gas-LDL critical point 
d for k B T/U A = 0.570, 0.580, 0.590, 0.595, 0.600, 0.610, 
0.620, 0.630. 

FIG. 2. The phase diagrams, with coexistence lines and 
critical points resulting from MD simulations. (a,b) P-T 
phase diagram. Panel (b) is a blow up of panel (a) in the 
vicinity of C\. Circles represent points on the coexistence 
fines: open circles are for the gas-LDL coexistence, filled cir- 
cles for the gas-HDL coexistence. Lines are guides for the 
eyes. The solid black line is the gas-HDL coexistence line. 
The red line is the gas-LDL coexistence line. The solid red 
line is stable, while the dashed red line is metastable, with 
respect to the HDL phase. The orange line is the LDL-HDL 
coexistence line. The triangle represents the triple point. The 
projection of the spinodal line is represented in (a) and (b) by 
diamonds with dashed lines. The spinodal line is folded in this 
projection, with two cusps corresponding to the two maxima 
in Fig.^. Critical points occur where the coexistence lines 
meet these cusps. The critical point C\ is for the gas-LDL 
transition, and Ci is for the gas-HDL transition. (c,d) P-p 
projections of panels (a,b). The colors, symbols and patterns 
of coexistence lines, triple and critical points are the same as 
in panel (a,b). Dashed blue, black and orange lines schemat- 
ically represent isotherms at the temperatures of Ci, Ci and 
of the triple point, respectively. Both gas-LDL and gas-HDL 
coexistence lines show a local maximum, representing the es- 
timates of C\ and Ci, respectively. In all the panels, where 
not shown, the errors are smaller than the symbol size. 
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